Anomalous scaling and Lee- Yang zeroes in Self-Organized Criticality. 



B. Cessac, J.L. Meunier,* 

We show that the generating functions of probabihty distributions in SOC models exhibit a 

T-H , Lee- Yang phenomenon Namely, their zeroes pinch the real axis at z = 1, as the system size 

' goes to infinity. This establish a new link between the classical theory of critical phenomena and 

' SOC. A scaling theory of the Lee- Yang zeroes is proposed in this setting. 

c^l ! 

^ ■ PACS number: 02.10.Jf, 02.904-p, 05.45.-f-b, 05.40.-fj 

o : 

(y-^ ■ In 1988, Bak, Tang and Wiesenfeld |Q] proposed for the first time a mechanism in which a dynamical system reaches 

' ' "spontaneously" a stationary state with some features reminiscent of a critical state. More precisely, by its only internal 

(-H ^ reorganization in reaction to (stationary) external perturbations, a system organizes into a state with scale invariance 

^ , and power law statistics. This effect, called Self-Organized Criticality (SOC), was quite unexpected, since usually, 

^ r the critical state of a thermodynamic system needs a fine tuning of some control parameter (temperature, magnetic 

' field, etc..) which is at first sight absent from the model introduced by BTW and from the many variants proposed 

' later |^,^. Furthermore, the stationary regime corresponds to a non-equilibrium situation where the (stationary) 

, fiux of external perturbation is dissipated in the bulk or at the boundaries, generating a constant flux through the 

S ■ system. As a consequence, one generally believes that the usual equilibrium statistical mechanics treatments using 



' ^ ' the concepts of Gibbs measure, free energy, etc ... cannot be applied for the analysis of SOC systems. 

c ■ 

O ' On the other hand, it is also believed that concepts like universality classes, critical exponents, order parameter, etc 



X 



borrowed from the equilibrium statistical mechanics of phase transitions are still relevant in SOC. Actually, the 



. identification of universality classes is one of the main goal in the SOC literature. However, since these concepts are 

> 

f — not defined via a thermodynamic analysis, alternative definitions are used. The dynamics of SOC systems occurring in 
chain reaction or "avalanche" like events, a set of observables characterizing the avalanches, size (s), duration (t), area 
[ (a), etc . . . are defined. Fix such an observable, say N, and compute the related probability Phin) Prob{N — n) at 
■ stationarity for a system of characteristic size L. The numerical simulations show that the graph of Pl{ji) exhibits a 
power law behavior over a finite range, with a cut-off corresponding to finite size effects. As L increases the power law 
range increases. This leads to the conjecture that as L ^ cx), Phin) converges to a probability distribution P*(n), with 
^ ' a power law tail having an exponent t„ called the critical exponent of the observable n. It seems commonly admitted 
C in the SOC community that a classification of the models can be made by the knowledge of their critical exponents 
O • ("universality classes"). Consequently, a large effort has been devoted to the computation of these quantities. 



Considerably less efforts have been made to establish a clear foundation of the basic SOC concepts and to clarify their 
connection to their classical statistical mechanics counterpart [^j5| . Clearly, this is a hard task since even preliminary 
notions like "state" or "thermodynamic limit", though intuitively clear, suffer from a lack of precise mathematical 
definition in this setting. In discrete automata like the BTW model ||] the state is the unique ergodic probability 
measure, /u^, of a discrete Markov chain, finite when L is finite. In continuous dynamical systems like the Zhang 
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model 0] there exists typically infinitely many ergodic measures and therefore one has to add additional constraints 
to define the state m a non ambiguous way [p|-p^. The probability distributions P^, for the observables a, s, t . . . 
are directly obtained from /i^ [slJlOl] but they contain less information. The observable a,s,t are simply indicators 
of the dynamics. There is no a priori reason to believe that the knowledge of Pl{s), PL{a), PL{t) or, even, of the 
joint probability PL{a,s,t) gives all the relevant (that is, allowing to classify the models into universality classes) 
informations about the state ^l- 

The thermodynamic limit L oo and the supposed "convergence" of /i^ to a "critical state" poses deeper problems 
since even the proof that there exists indeed a limiting state and that the probability distribution of avalanche 
observables are still defined in this limit remains to be done, for most models. The usual classical statistical mechanics 
constructions of the thermodynamic limit like the Dobrushin-Landford-Ruelle's (DLR) cannot be directly applied 
because of the a priori absence of a Gibbs formalism. On the other hand, the methods used in interacting particle 
systems, allowing to define the dynamics in the thermodynamic limit, on the basis of the Hille-Yoshida theorem and 
the properties of Feller processes, requires locality conditions which are broken in SOC model. One has then to 
develop new ideas for non-Fellerian Markov processes and this has only been done in a few examples [l^ . 

But one of the main problem is the treatment of the data obtained from finite size systems simulations itself and the 
extrapolation to the L oo limit. Indeed, though it was believed in the earlier SOC papers that this extrapolation 
can be handled by classical finite-size scaling |l3), further investigations proposed alternative scaling Jl^-p^ and, at 
the moment, there is no agreement on which scaling form applies. Consequently, a lot of efforts have yet to be devoted 
to the understanding and analysis of SOC models. 



Though the analogy between self-organized criticality and usual critical phenomena is the core of the SOC paradigm, 
it is remarkable that, up to now, some well developed techniques of analysis of critical phenomena have not been 
adapted to the study of SOC models. A phase transition has different manifestation. It is in particular characterized 
by a singularity of the thermodynamic potential (free energy, pressure). At a phase transition point, and for suitable 
interactions, the free energy, which is the generating function of the cumulants, exists in the thermodynamic limit but 
it is not analytic (in a fc th order phase transition it is C*^^^ but not C*^). In many examples, the failure of analyticity 
is manifested by the Lee- Yang phenomenon IT^ ]. For finite size systems the partition function is a polynomial in 
a variable z which typically depends on control parameters like the temperature or the external field. Since all 
coefficients are positive there is no zero on the positive real axis. However, in the thermodynamic limit, at the critical 
point, some zeroes pinch the real axis at z = 1, leading to a singularity in the free energy. The finite-size scaling 
properties of the leading zeroes and of the density of zeroes near to z = 1 determine the order of the transition 
and also the critical exponents in the case of a second order phase transition ||l9| ] . 

A natural question is whether there exists a similar property in SOC, namely can we exhibit a' "free-energy like" 
function, developing singularities in a similar way in the infinite lattice size limit. Though there exists a huge literature 
about the Lee- Yang zeroes, there is, to the best of our knowledge, no attempt to study Self-Organized Criticality from 
this point of view. In this paper, we show that the cumulants generating function of the probability distribution of 
the observables a, s, t, . . . have this property. More precisely, the expected convergence of Pl to a power law induces 
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a Lee- Yang phenomenon for the corresponding cumulants generating function (eq. |l|). We show that this effect is 
related to the observed divergence of the moments. Furthermore, a scahng theory of the Lee- Yang zeroes is proposed. 

After some prehminaries (section , we give exphcit analytical results (section ||) in several cases used as guidelines 
for subsequent analysis of a SOC model (section We first study the truncated power law case where the cut- 
off tends to infinity when a parameter L, (corresponding the lattice size in SOC models), tends to infinity (section 



II A). We give in particular an analytic expression for the zeroes. Then, we investigate the effect of a smooth cut off 



(section II B). We first discuss the properties that this cut off must have, extrapolated from numerical simulations, 



and present some of the ansatz found in the literature (section II B 1 ) . We then explicitly compute the Lec-Yang 
zeroes for a probability distribution obeying the finite-size scaling form proposed in |ll3|] and converging to a power 
law as i — > cxD (section II B 2| ) . We show in particular that, when the power law exponent r is larger that 1 there is 



a violation of the scaling usually observed in classical critical phenomena; namely, there is an anomalous logarithmic 
dependence on L for the angle that the zeroes do with the real axis in the t = log(z) plane. We also show that 
when T > 1 a bias is artificially induced by the numerical simulations, when the size of the sample used to generate 
the empirical probability distribution is fixed independently of the lattice size. This effect can be analyzed with the 
Lee- Yang zeroes (section IIB 3| ). We then briefly study some other scaling form proposed in the literature and the 



effect of a finite size scaling violation on the Lee- Yang zeroes (section IIB 4). Finally, in the last section, we present 
numerical simulations for the Lee- Yang zeroes in the Zhang SOC model and compare them to the theoretical results 
obtained in section We see no clear cut evidence of finite size scaling violation, but show that this model is quite 
sensitive to the numerical cut off induced by a lattice size independent sampling. This can raise some doubts on the 
conclusions about scaling (Finite Size Scaling, multifractal or whatsoever) which can be drawn from some large lattice 
simulations done on this model in the literature. 

I. PROBABILITY DISTRIBUTION AND LEE- YANG ZEROES. 
A. The finite size system. 

Let Phin) = Prob{N = n) be the probability distribution of the avalanche observable N G 1 . . where the index 
L refers to the characteristic size of the system, ^l, the maximal value that N takes is finite, whenever L < oo, but 
diverges as L oo. Therefore, the fmiction: 

ZL{z) = Y,z-PL{n) (1) 

n=l 

where z e C, is a polynomial of degree ^l. In particular, since Zl{z) is an analytic function of z in the complex plane, 
all its moments exist. Denote by E[]l the expectation with respect to PL{n). Call : 



mi(g) = ^Pi(,i)n'°i'i?[n«]i (2) 



dcf 
n=l 

where q is a real (positive) number. For integer the niL{qys are the moments of P^^n). Note that the normalization 
of PL{n) imposes Zl{1) = ^^(0) = 1. 



3 



For finite L, Zl{z) has zeroes in C, that are either real < 0, or complex conjugate. Denote them by ZL{k), k = 
1 . . .£,L and order them such that < — 1| < . . . < \zL{k) — 1| < ■ • • < — 1|- Note that 2 = is a trivial 

zero, of multiplicity one, since Pl(1) > 0. Write ^^(fc) = RL{k)e^^'^'--''^ = l + ri(fc)e*'^^('^'. Since all Phin) are positive, 
Zl{z) has no zero on the positive real axis for finite L. Consequently, the log-generating function^ log[ZL{z)] is well 
defined on IRl. Furthermore : 



GUt) ""^^ log[ZL{e')] (3) 



is an analytic function of t. Call 



xUq) = —log[ZL{z)] 



(4) 



where z = e*. The quantities xl{<i) are easily expressed in terms of the Lee- Yang zeroes : 



B. The "thermodynamic" limit L ^ 00. 

1. Divergence of the moments and Lee-Yang phenomenon. 

As already written in the introduction, a mathematical definition of the thermodynamic limit in SOC is a difficult 
task, beyond the scope of this paper. However, in [ pO| , we developed a dynamical system approach for the Zhang 
model. Then, the thermodynamic formalism can be used to define the finite size SOC state of as a Gibbs 

measure^ in this setting. It is then shown that the joint avalanche size distribution, for example, can be obtained 
in this formalism via a proper potential. The corresponding generating function for the time correlations, called the 
topological pressure, is the formal analog to the free energy. In this setting, it is argued that the critical behavior 
expected in the thermodynamic limit, is manifested by a non analyticity of the topological pressure as L ^ 00, that 
can be linked to the loss of hypcrbolicity characterizing the limit L ^ 00 of the Zhang model The loss of 

analyticity can be easily detected by looking at the generating function (^. Indeed, its zeroes exhibit a Lee- Yang 
phenomenon. 

The paper |20| is devoted to dynamical system aspects and to the mathematical foundation of a thermodynamic 
formalism for the Zhang SOC model, and the link between the scaling theory of Lee- Yang zeroes in classical critical 
phenomena and general SOC model is not addressed. This is the aim of the present paper. The results developed 
here are therefore complementary to pO| but are independent. 



^There is obviously a formal analogy between (^) (resp. (^) and a partition function (resp. a free energy). 

^The particular structure of the Zhang model allows to symbolically encode the dynamics. In the framework of the thermo- 
dynamic formalism a Gibbs measure is a probability measure weighting the symbolic chains encoding the trajectories with an 
exponential weight called a potential (see liB] for details). 
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The present paper focus on the analytic properties of the fog generating functions of the probabiHty of avalanches 
indicators, when L oo. It intends to analyze the variations in the Lee- Yang zeroes properties if one uses the different 
scaling forms found in the SOC literature. Consequently, we collected the minimal implicit assumptions used in the 
SOC literature and we infer the consequences they lead to. This means that the result developed a priori hold for all 
SOC models. 

It is first assumed that PL(n) converges to some probability distribution P*{n), n = 1 . . . oo. It is furthermore 
assumed that P*in) has a power law tail^, namely P*{n) ~ for a certain n range, n = no . . . oo where ng < oo, VL. 
The number no depends on the model and on the observable and introduces an extra parameter in the characteristics 
of the probability distribution. In the computations done in this paper there is no loss of generality is assuming that 
no = 1. Therefore, in the sequel, P*[n) will stand for the power law n = 1 . . . oo. 

The measured exponent t in SOC belongs to the interval ]1,2[. K — P*{1) is the normalization constant. Conse- 
quently, K = where C is the Riemann C, function Under the above assumptions, the moments mL{q) behaves 
asymptotically like X^nJ^ n^^"^. This sum diverges for all g > t — 1. It is numerically observed that mL{q) diverges 
like mL{q) ^ L'^^'^K A central issue is to compute the scaling exponents given by: 

, del \og{mL{q)) log(xL((7)) , . 

a(q) = hm — - — — — = hm — — — — (6) 

L^oo log(i) L^oo log(L) 

(j{q) is an non decreasing function. Its Legendre transform is found under the name of "multifractal spectrum" in the 
SOC literature p^ ] though it has no direct connexion with the fractal geometry of the invariant set. 
Since P* (n) is a probability distribution the limiting generating function: 

oo 

Z*(z) = lim = y P*(n)z" (7) 

11=1 

is still an analytic function in the open unit disc in C. However, the log-generating function of P* (n) is not analytic near 
to z = 1 since the derivative of order q > t — 1 > diverge. The corresponding singularity is related to the behavior 
of the zeroes in the vicinity of z = 1. More precisely, fix e > arbitrary small, call /^(e) — {i \ {zLii) ~ 1| < e} and 
ni(e) = #/i(e) where # denotes the cardinality of a set. Then the divergence of is governed by the zeroes 

which accumulate in Namely, the sum (|^) contains a singular term : 

,.(L,.,,).(-ir^2.(,-l)!^S^^5i^ (8) 

k=l I 

which diverges as L ^ oo, while the remaining part in the sum is regular and is bounded by '^'^ J^^' as L ~> oo. 



^Note that the limiting probability distribution is defined only if r > 1. 
*In general the normalization constant depends on no. 
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2. Scaling of the zeroes in classical critical phenomena. 



In the theory of classical critical phenomena, it is possible to relate the scaling exponents of quantities such as 
magnetization or latent heat, susceptibility, etc ... to the behavior of the Lee- Yang zeroes near to z = 1. There exists 
a scaling theory based on earlier works by Lee and Yang Grossmann and Rosenhauer ||2^, Abe Suzuki 
[psf , Privman and Fisher |2^], Itzykson et al. ||3^, Glasser et al. Many analytical and rigorous results are also 
known (for example [^2[-p5|). A lot of efforts have been devoted to the study of ferromagnetic systems (e.g. Ising or 
Potts models) though many other examples have also been studied in the literature. In this setting, one distinguishes 
the zeroes in the complex magnetic field (called Lee- Yang zeroes) from the zeroes in the temperature plane (Fisher 
zeroes). In the first case, the zeroes lie on the unit circle for a large class of models including the Ising's one. 

The Fisher zeroes usually approach t = Q in the t — log(z)-complex plane with a constant angle <j) (this is the case 
for the Ising model and mean- field ferromagnetic models ) . This allows to obtain simple scaling expression for the 
singular part of the free energy ^{t) where t is the reduced temperature. In this setting, an analytic expression for 
/*(t) has been obtained by Grossmann and Rosenhauer and, later on, extended by Itzykson et al. js^ by using 
the renormalization group theory. This approach has been extended by Glasser et al. |]3l]] to mean-field models. In 
the thermodynamic limit ff- ^ A±\t\^^°' where A± are universal constants (± label the two magnetic phases at low 
temperature), and a is the critical exponent for the specific heat. It follows from the renormalization group analysis 
p9| that the singular part of the free energy obeys a Finite-Size Scaling form: 

rit,v) = ^mAv)^] (9) 

where V is the finite volume and J- a universal function. Accordingly, the n first Fisher zeroes are given by : 

27r 1 



tv{n) = 



{Al +Al- A+A^ cos(7ra)) V 



cos(7ra)— — 

The angle (p is related to A±,a by tan [(2 — a)(f)] — - — gi^^^^) ^ ■ This situation, where the zeroes approach the 
singularity with a constant angle 4> and where the modulus scales like the volume to a certain power will be referred 
to as normal scaling in the sequel. 

It seems a general observation p6| that the zeroes lie on a curve or a union of curves dividing the complex plane 
in different regions of analyticity of Z*{z), corresponding to different phases. More precisely, it has been recently 
proved by Biskup et al. |Q] that the zeroes lie on curves with a simple analytic expression and accumulate in the 
thermodynamic limit on loci where the various branches of the free energy have the same modulus. This last result 
suggests that a wide extension of the Lee- Yang phenomena can be made toward dynamical systems near to a critical 
point. We now develop this aspect for the analysis of the log generating function of probability distribution in the 
SOC framework. In the sequel, we will not distinguish between the Lee- Yang zeroes and the Fisher zeroes and we 
will use the generic terminology "Lee- Yang" for the zeroes. 
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II. SCALING THEORY OF SOC AND LEE- YANG ZEROES. 



In this section, we establish analytical results for various finite size scaling forms found in the SOC literature. These 



results are then used in section [II for the analysis of the empirical data obtained from a numerical simulation of a 
SOC model. As a matter of fact, for finite size SOC systems, the power law is truncated by a cut off characterized 
by a length scale A^, usually different from ^l. The starting point is therefore the analysis of a truncated power law 
with a sharp cut-off at a value A^. This is a useful example for subsequent analysis since the analytic form of PL{n) 
and the cut-off is known. 

A. Zeroes of a truncated power law. 

Assume that Phin) — — ^■■■S.l where Cl is a normalization constant and = A^ — L^,t > l,/3 > 0. 

Furthermore, assume that r < 2. 

1. Scaling of the moments and log generating function. 

For < g < T — 1, rriL^q) '^''ciT)'^ consequently cr{q) — 0. The non-zero scaling exponents cr(q) can be 
obtained from the following integral approximation of mL{q), which becomes exact in the limit L ^ oo, provided 
q > T - 1 : 

Cl 
1 

Then, a{q) — (3{q + 1 — t) for (real) g > t — 1. Note however, that for finite size, one has additional L dependent 
terms which have to be considered when extrapolating from numerical simulations. It is also interesting to note that 
formula ( pj]) gives useful informations on the rate of convergence of mL{q) to a constant for g < r — 1. Indeed, the 
convergence is not uniform in namely the closer is g to r — 1 the slower is the convergence rate. This means 
that a systematic bias due to finite size is introduced in the numerical simulations when dealing with the g's close 
to T — 1. This produces a spurious curvature, near to r — 1, for the function cr{q) extrapolated from numerical data. 
This effect, which disappears as L — > oo, can lead to misleading conclusion since it can be interpreted as an evidence 
of a multifractal scaling. 

The scaling exponents can also be obtained from the scaling of the log generating function. Indeed, the generating 
function writes: 



niLiq) ^ CiAi+i-^ / u'^-^du = ( Afi-^ - 1) (11) 



^lW = 1 + E PLin)ie'" - 1) ^ 1 + CLA'r [[ ^'"( 

1 T 



Set t' ^ K^t and 



•"^ n=l ' n=l 



m / .-(e*- - l)du = > : - / u-^du = ^ (13) 
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Note that since r < 2, this integral is finite as can easily be checked by integration by part. Therefore the 
commutation of the integral and the series is allowed. 
Consequently, 

ZLit) ^ 1 + CLAl-^yjit') - 1 + CLL^'^^-^mtL^) (14) 

and : 

GL{t) = log(l + CLLf"^^--H{tLf')) (15) 

ip{t) is a smooth function of t which vanishes as i ^ 0. Therefore, for t : 

GL{t) ^ CLL^^'^^-^mtL^) (16) 

which gives the right scaling for the moments by differentiating with respect to t at < = 0. One remarks that this 
scaling form is analogous to the form (^. 

2. Lee-Yanq zeroes. 



From equation ( p[ ) the zeroes are approximately given by: 

i^it') = -ClKI-^ = -ClL"^^-'^ (17) 

Since r > 1, A]^~^ diverges, ipit') is an increasing function of the real variable t' which vanishes as t' = and tends 
to —oo when t' — oo. Furthermore, for any K > 0, tp{t') is bounded by V'(-f^) in the ball \t'\ < K in the complex 
plane. Consequently, the equation ( pT| ) can be fulfilled only if t^(fc)'s have a diverging modulus as L grows. On the 
other hand, since tL{k) — converges to zero, \t'j^{k)\ must grow slower than A^. It grows in fact like log(Ai) as 
shown below. 

Note that, conversely, when r < 1, —ClL^^'^~^^ goes to zero in the thermodynamic limit |. Then, the zeroes are 
formally given by : 



where ijj'^^ is the k th branch of the inverse of ijj in the complex plane. Consequently, the zeroes have a simple scaling 
in these case similar to (p^). In particular, the argument does not depend on L. 

The observed values of the critical exponents in SOC, t > 1, induces anomalous scaling that can be observed on 
the Lee- Yang zeroes. Though the form dlj) can be used to compute the Lee- Yang zeroes, it is easier to use: 



^Though the limit probability is not defined it is nevertheless possible to investigate the finite size scaling properties for the 
zeroes of the finite size generating function. 
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Zdt) = J2 ^iWe*" ClAI-^ / ^ h{u,t' 



)du 



(18) 



where t') = ' 



There exists several techniques to compute the Lee Yang zeroes in statistical mechanics. A 



standard way is to argue that the asymptotic free energy admits different analytic continuation in different regions of 
the complex plane, separated by Stokes lines where the zero accumulate in the thermodynamic limit. Indeed, because 
of the large number of terms in the polynomial which make up the partition function, the behavior tends to become 
dominated by some set of the coefficients. Thus we have different analytic functions in different regions of the complex 
plane. These functions have oscillating phases but smoothly varying amplitude. The zeroes locates then on Stokes 
boundaries where two types of behavior have comparable magnitude ^0 31 3^. The Stokes boundaries becomes cuts 
in the thermodynamic limit. Across the boundaries the free energy has a regular real part and jumps in the imaginary 
part p^. 

Applying this strategy to our formal partition function ([T^), one identifies easily two regions. For real t' , as u 
grows from to to oo, h{u,t') first decay like u^'^ until a minimum m_ = p-. Therefore, u_ > ^ when t < t. For 
u > h{u,t') grows exponentially like e* Therefore, when t' is small the integral in (^8|) is essentially dominated 
by the algebraic decay u^'^ and f\ h{u,t')du -^^^ \A^i7^ ^ l] ■ On the other hand, for large i', u_ ^ and the 
algebraic part is negligible compared to the exponential part. Hence / i h{u,t')du ~ e* "du — e* — e'^ . This 
argumentation extends to complex t' and suggests that one can roughly divide the t complex plane into two regions 
where ZL{t) has a different analytic form: for sufficiently small t' the algebraic part dominates, while for large t' the 
exponential part is dominant. Then the zeroes have to stay at the place where the two forms are of the same order. 
Therefore an approximate equation for the location of the zeroes is given by: 

1 r 



t' 



3*' + ait' - a-i 



ATI 

T-1 



(19) 



where a\ 



l-T 



02 



The solutions of ^ = 



a:- 



— are given by : 

iL(fc) = Aiti(fc) 



(20) 



where Wk{x) is the k th branch of the Lambert function [p6[ . Note that the Lambert function has infinitely many 

t' A^-i 

branch and consequently the equation ^ — — has infinitely many solutions. Indeed, in replacing the initial sum 
by an integral in ( |l8|) we have introduced spurious zeroes which have to be removed for finite L. In the sum ( [l^ ) one 
has a step of integration which defines a lower cut-off in the scales one has to consider. Consequently, only the 
branches fc = — ^ . . . -j, where one takes into account the symmetry of the zeroes with respect to the real axis, are 
relevant. 

The Lambert function W-k admits the following expansion, for large | log(x)| pGf : 

(loglogfe(a;))" 



Wk{x) = logfc(a:) - log(logfc(a:)) + ^ ^ 



(>0 r?j.>l 



(logfe ^) 



(21) 
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where logj, is the k th branch of the complex logarithm and Qm = [1+™] is expressed in terms of the Stirling 

numbers of the first kind (-1)™+" [;^] [g. 

The double series J2i>o J2m>i ^\°iog°.%^'+" absolutely convergent for sufficiently large | log(a;)| [Q. Since we 
are only interested in the asymptotic divergence when L grows, one can therefore neglect the series in the asymptotic. 
Note however that the convergence to the asymptotic regime where the series becomes negligible is faster when the 
product f3{T — 1) is larger. 

The term log(log^(^^^— j-)) cannot be neglected compared to ^ogki-^tr) since it contains crucial k dependence for 
the real part of tL{k) (see below). It is interesting to note that it introduces a log(logi)) finite size scaling correction. 
A similar correction as been found in |^ for the Potts model with g > 4. 

The corrections due to the other terms of eq. (|l^) become rapidly negligible as can easily be seen by a perturbation 
expansion. 



One finally obtains the following asymptotic form for the Lee Yang zeroes of the truncated power-law: 

log(f^) + i log (\0g\^) + 4fc2^2^ 

Q(ti(fc)) 



2kTT 

a7 



■ arctan 



Al 
27rfc 



(22) 
(23) 



The term log(logfc( ^^Ji )) introduces a k dependence which implies in particular that the zeroes (in the z plane) 
do not lie on circle but on a more complicated curve (see Fig. This dependence remains important, even for the 
first zeroes, up, to very large L, especially if r is close to 1. Indeed, for a fixed k the term log^ ( ) dominates the 



term 4fc^7r^ only for L >> ((r — 1) 



(say, for r = 1.25, /3 — 2.67 this corresponds to a L >> 1500). 



The arctan term in the imaginary part acts essentially as a phase term — arctan 



which is slowly 



2/C7I 



Furthermore, since the arctan is bounded above by 



varying (in the k variable) compared to the dominant term 
■| it is rapidly negligible as k grows. Therefore, one can consider with a good approximation that 3(iL(fc)) ^ 

The argument of t^lk) formally corresponds to the angle that the Fisher zeroes do with the real axis in critical 
phenomena. For a power law with r < 1 this angle is independent of L as discussed above. Conversely, for r > 1 it 
is given by : 



/ 



Arg{tL{k)) ^ arctan 



2fc7r 



log(^) + i log (log'(f^) + 4fc2^2) 



(24) 



One observes therefore a logarithmic deviation to the normal scaling on the argument of the tL{k) 's. The conclusion 
is therefore that, though the truncated power law obeys the classical finite size scaling form (^), the Lee- Yang zeroes 
display nevertheless an anomalous scaling due to the exponent r > 1 (see Fig. ||). 

In the z plane the zeroes ZL{k) = e*^'-'^^ are approximately given by: 



RUk) = |zL(fc)| - l + 5R(tL(fc)) 

„ /,N 2fc7r 2fc7r 

0L{k) = Q(iL(fc)) ^ — = — 



(25) 
(26) 
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Therefore the arguments O^^k) of the zeroes in the z complex plane are uniformly distributed in [— tt, tt] with a 
good approximation. 

Finally, one can determine the exponents r, f3 from the Lee- Yang zeroes. The exponent f3 corresponds to the scaling 
exponent of the correlation length ^l. The equation ( p6| ) provides a straightforward way to compute it. Furthermore, 
the exponent r can be obtained from eq. (p2|). The term 47r^fc^ certainly rapidly dominates the term log^ ( ) in 
the modulus i?L(fc) as k grows. This is a fortiori true for k Al which correspond to the zeroes the farthest from 
z — 1. Consequently: 



RLiO) 



r- 1 



1 




1 










r- 1_ 





If one takes k 



Ai. 

2 



(corresponding to an angle tt), one has: 

Al log(i?LW) 



lim 

L — *oo 



(27) 



(28) 



log(Ai) 

which allows a possible determination of t from the scaling of the Lee- Yang zeroes. Note however that the convergence 
is logarithmic in L. 



3. Numerical checks. 



Since it is easy to generate numerically a power law distribution one is a priori free to choose any values for r and 
f3. However, the closer r is to 1 the slower is the convergence to the asymptotic regime where the formulas obtained 
in the previous section hold. More precisely, the rate of convergence is essentially governed by the product (3{t — 1). 
The closer is r to 1 the larger has to be f3. But the larger /3, the faster the degree of the polynomial increases with 
L and therefore the time needed for the computation of the zeroes increases. On the other hand, since the theory 
developed here is independent of the actual value of /3,t (provided r > 1), we mainly studied examples where [3 = 2 
and (3 — 2.2 which gives a reasonable increase in the polynom degree, and t ~ 1.9 such that the product f3{T — 1) ^ 2. 

We first depicted the pattern of zeroes, for different sizes, in Fig. |l|. 




near to 
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One notices the slow convergence to the unit circle, and the shape of the curve near to z = 1: this is not a circle. 

We plotted Fig. || the real and imaginary part of the zeroes in the t plane. For the real part we tried a fit of the form 
r(fc) = - — '°s(°^+2 i°g|i°s (a)+47r fc ]) ^j^^^^ a is a free parameter. The formula ( [2^ ) gives ath — but, in the several 
approximations we made, we neglected some constants, and one expects a to be different from the theoretical value, 
with an error that should decrease as L grows. The result of the fitting is represented Fig. ^ for the 200 first zeroes. 
We found indeed that the experimental value aexp is closer and closer to its theoretical value as L increases and that 
our approximation is better and better as L increases. For L — 50, ath = 7.87 x 10"'', aexp = 5.6 x 10^^ ±4 x 10^^; for 



L = 100, ath = 2.26 X 10" 



1 ^''exp 



= 1.5 X 10"''±3 X IQ-^ and for L = 150, ath = 1-09 x 10" 



1 ^''exp 



= 7.1xl0-5±810 



-7 



The imaginary part is plotted fig. § b together with the theoretical prediction 3(ti(A;)) = (straight fines). We 
remarked a slight deviation of the first zeroes to ^ due to the correction term in (p6|). 
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FIG. 2. Lee- Yang zeroes in the t complex plane, versus k for various L values. Fig. Real part multiplied by Al. Fig. ^ 
b: Imaginary part. The fitting curves are plotted in color, full lines. 



Fig. ^ a, we plotted the argument of t^^k) as a function of L for fc = 1 . . . 5. We notice the logarithmic de- 



viation to the normal scaling as predicted by formula (24). We tried to fit these curves with a fit of the form 



arctan 



2kn 



where a, 7 are free parameters. In fig. ^jb, we show the different scaling 



occurring for r < 1 (normal scaling) or r > 1 (anomalous scaling). 
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FIG. 3. a Argument of the Lee- Yang zeroes in the t complex plane, versus L for various k values, P — 2,t — 1.9. Fig. ^b. 
Normal (r = 0.2) and anomalous (r = 1.21, t = 1.9) scaling of the angle with the real axis for various r value. The fitting 
curves are plotted in color, full lines 



Finally we argued above that the exponents f3 and r can be determined with a good accuracy from the scaling of the 
zeroes. In fig. a we plotted Arg{zL{5)) — 6l{5) as function of L ior (3 — 1.5, r — 1.6; /3 — 2,t— 1.21; /3 — 2.2, t — 1.9 
and L = 20 . . . 120. We choose the 5th zero rather than the first ones since for the first zeroes the correction coming 
from the arctan term in eq. (p3|) influences slightly the scaling for small L. We tried a fit of the form where 



Cth ^ lOn = 31.415. 
2.205 ± xO.OOl 



The fits shown on the figure gave respectively : (3 = 1.507 ± 0.002; /3 = 2.007 ± 9 x 10"^; [3 ■ 



For the determination of r we used eq. ( |27[ ) for r = 0.2; 1.21; 1.9 and /? = 2. We have interpolated the data 
with a fit form e/5'^i°g(''^)/^^ where a,T were free parameters. For L = 20 . . . 120, we found respectively 0.1998(2) ± 
6.10"^, 1.207(5) ± 7.10~^, 1.89(5) ± 0.0001 for the value of t. This is satisfactory especially when taking into account 
the smallness of the L's we considered. 
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FIG. 4. a Argument of the fifth Lee- Yang zero in the z complex plane, versus L for various [3 values. Fig. Rl{t^) as a 
function of L for various r values. The fitting curves are plotted in color full lines. 



B. General case : the structure of the cut off. 

The observed probability distribution of avalanches observables in SOC models is in general a power-law truncated 
by a cut-off function associated to finite-size effects. Except for a few case |3^], the analytic form of the cut-off is 
not known. Consequently, various scaling form have been proposed in the SOC literature. In this section we grasp 
the most common scaling forms and discuss their effect on the Lee- Yang zeroes. We also investigate the effect of the 
sampling cut-off inherent to numerical simulations. 

1. General assumptions about the cut off function. 

Since Phin) is expected to converge to a power law Kn~'^ one can write, without loss of generality, the finite size 
probability under the form; 
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PlH = ^fUn), 1 < n < a (29) 

where Cl is a normahzation constant depending on L. fLin) is the finite size cut-off. 

The graph obtained from numerical simulations suggests that fhin) is a regular function which obeys: 

(i) lim fL{n) = 1, Vn<a- 

L — yoo 

(a) Vp > 0, lim nPfLin) ^ 0, VL 

n — >oc 

The property (i) corresponds to the pointwise convergence to a power law. (ii) characterizes the observed fact that 
the tail of PL{n) decreases faster than any power of n (e.g. is least exponentially decreasing). 

These properties are not sufficient for a scaling theory and further assumptions have to be made. We now discuss 
various scaling form that one can find in the literature but also the numerically induced cut-off effects. 

Finite-Size Scaling. In 1990, Kadanoff et al. jl^ proposed a Finite-Size Scaling Ansatz where : 

fL{n)=g{—) (30) 
g being a universal function. = is the characteristic scale for a lattice of size L and = aA^ where a is a 



constant. The case of a truncated power law developed in section II A corresponds to the particular case where a = 1, 
and g{x) is equal to 1 for x G [0, 1] and zero otherwise. Note that in this case the property (ii) has the consequence 
that: 

{iii) lim g{x) = 1 

multifractal scaling. In the same paper, Kadanoff et al. |l^ discussed another form of scaling, that is: 

\og{PL{n)) _J \og{^j \ 

where / is universal and does not depend explicitly on L. Lq, uq are some constants that can be omitted by a suitable 
redefinition of the quantities, then: 

PL(n) =L^('°s"/'°9-^) (32) 

log(-!!-) 

This representation is called by the authors an / — a representation, a being the quantity j^^^ j? ^ . In this case one 
has a whole spectrum of scaling indexes, i.e. all the values taken on by In the general case / is non linear. Then 
the universality class is given by the function /, rather than by a finite set of critical exponents. In this case, the 
scaling exponents are non linear functions of q. 

Finite Size Scaling violation and convergence to a power law. Another scaling form has been introduced by Lise 
and Paczuski [ p^ from numerical simulations on the OFC model |^. Lise and Paczuski analyzed their data with the 
following form for Pi(n): 



B.C. is very grateful to M. Paczuski for illuminating discussions on this topic in Bielefeld. 
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P^(n) = CLn-^L^''(^^); n = 1.. .L^^ (33) 

where Pl is now L dependent, [3l < P < oo and (3l —>■ (3 as L ^ oo. Furthermore, the numerical plot of Fl{x) 
in |l6j suggests that Fl(x) converges to a "step" function Y{x — /3) as L — » oo, where Y{u) = 0, u g] — oo, 0[, 
y(0) — C and Y{u) = — oo, u s]0, oo]. The finite size scaling case corresponds to Fl{x) — '"^^^^g^ — ^ where g is 
defined in (|3^). This example is quite interesting since it gives an example of a probability distribution violating ( [30| ) 
but converging nevertheless to a power law (namely the exponents r and (3 are still meaningful). We remark indeed 
that the corresponding probability distribution is not multifractal in the sense of since it is easy to check that 
the scaling exponents are given by a{q) — (3{q + 1 — r), Vq > r — 1 and (j{q) — 0,\/q < t — 1. This case is therefore 
intermediate between the Finite-Size Scaling (pCl) and the multifractal case (pi). 



In the following sections we discuss the behavior of the Lee- Yang zeroes in these different cases. 

2. Finite-Size Scaling. 

We show in this section that the behavior of the Lee- Yang zeroes in the finite size scaling case is essentially the 
same as for the truncated power law, provided g in ( ^0| ) fulfills the conditions (i),(ii),(iii) in the previous section. 

It is well known that the moments obey the same scaling. This can be recovered from the generating function (|^). 
Provided r < 2 it scales as i ^ oo like : 

ZL{t) ^ 1 + ClAI-^ ho^it') - T 1 (t')l - 1 + CLAi-^T„(t') + ^ 1 + CLA[-^T„(i') (34) 

for t — !■ in the t complex plane. In this equation t' — A^t and : 



T,{t') / u-^g{u){e* " - l)du = J] ^ / u^-^g{u)dv 
Jo Jo 



(35) 



Consequently, near to t = 0, ZL(t) - 1 + CLL'^^^^^^T^itL'^) and 6^(0 - CLL^^^-'^^Tc^itLP). As expected one 
obtains the same scaling as for a truncated power law. The only change is that the function T depends now on the g 
function. 

The zeroes of Zl [t] are therefore well approximated by : 

Ta(t') - -ClA^-1 (36) 



which is completely analogous to (17). Therefore, the same conclusions hold: the equation (|36|) can be fulfilled only 



if the solutions, i^(fc), have a diverging modulus as L grows. The computation of the zeroes is essentially the same 



as in section (II A) which slight complications due to the presence of the function g. 



One can write : 

ZL{t)'^CLA\;^ I h{u,t')du (37) 
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where now h{u, t') = u~'^ g{u)e* Recall that, by hypothesis, g{u) is an at least exponentially decreasing function. As 
u grows from to to oo, h{u,t') first decay like u^'^ until a minimum u_ after which h{u,t') grows exponentially like 
U- tends to as i' — s- oo. More precisely, if t' is sufficiently large, from (iii) g{u) is essentially 1 on the interval 



u 



[0, M_] and therefore u_ is approximately given by t'uZ^ = tuZ'^^ ■ Consequently u_ ^ p-. When u is sufficiently 
large the decay coming from g{u) compensates the exponential increase of e* Therefore, there is a maximum 
after which h{u,t') tends to zero with a rate given by g{u). Here, h{u,t') is essentially (7(u)e* Therefore u+ is 
given by t'g{u^) + g'{u+) — or t' ~ — , where ~ is a function which increases faster than u by 

(m). Consequently, u+ diverges as t' — > oo. Since a is bounded by assumption, u+ > a for t' (resp. L) sufficiently 
large. But, from eq. (36) we know that the modulus of the t' corresponding to the zeroes diverge as L — > oo. 
Consequently, provided L for is sufficiently large, the zeroes lie in a region where it+ > a. Guided by the wisdom 



coming from section II A we can assume that the zeroes accumulate onto a curve 7l which separates the complex 
plane into two regions. In the first one, the behavior is dominated by the algebraic part and the integral in ( ^ ) is 
essentially f" h{u,t')du ^ [-'^L^^ ~ a'^^^]. In the second region, the exponential part dominates. Provided t' is 
sufRciently large >> a), the variations of u^"^ g{u) are small compared to the increase of e* " and this function 
can be approximated by some constant Tg. Hence J" h(u^t')du ^ Tg J" e* "c?w 
Consequently, for sufficiently large L, the zeroes are well approximated by 

1 



e''" + ait' 



02 



(38) 



where ai — p" ,02 e ' . 

This equation is similar to eq. (^). The zeroes are now given by 



tL{k) 



aAi 



(39) 



Consequently, one finds that the pattern of zeroes in the finite-size scaling case is essentially the same as the power 
law case, up to a correction depending on a,Tg. More precisely, using the series expansion ( ^l| ) of the Lambert 
function one finds : 



log( 



£l\(r-l) 



ilOg (lOg^(- 



aAr 



2kTi 
aAr 



(40) 
(41) 



where k = 1 . . . ~ aA^. 

In the z plane this essentially results in a trivial re scaling of the argument 3(ti(fc)) and a slight change in the 
modulus i?L(fc) : 



(42) 



where the superscript FSS (resp. PL) refers to the Finite Size scaling (resp. Power Law) situation, kl ^ ("T^) ""^^ 
is therefore a scaling factor which tends to 1 as expected since the zeroes have to accumulate on the unit circle. In 
equation (E2|) we neglected the more complicated dependence coming from the log(log) term. This has a small effect 
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on the first zeroes but becomes negligible as k grows. This provides a way to determine kl. Setting R£^{tt) for the 
farthest zero from z = I: 



The argument of ti(fc)'s is 



(43) 



Arg{tL{k)) ^ arctan 



2fc7r 



[- \og{^^^j^) + \ log (log2(fiI^) + Ak^^^) 
The cut off g modifies therefore the value of the angles but not the scaling. 



(44) 



We numerically checked these result in the following case. We generated a probability distribution given by 

PL{n) = CLn-^g{^)- n^l...aL^ 

where : 



(45) 



(46) 



We fixed the parameters to the values: /3 = 2, t = 1.9, 7 ~ 2.5, a — I. We show Fig. ^ the collapse of the curve of 
zeroes to the corresponding curve for a power law with the same r, /3. k was computed from the ratio (^3|). We found 
K = 1.0001 for L = 100. 



0.04 



0.02 



■0.02 



■0.04 



POW'k 

"""«».. FSS 




1 1.0002 1.0004 1.0006 1.0008 1.001 1.0012 1.0014 
Re(z) 



FIG. 5. Pattern of zeroes in the z complex plane for the power law case, l3 — 2,t = 1.9, L — 100, where the real and 
imaginary part have been multiplied by k, and for the FSS case P — 2,t = 1.9, 7 = 2.5, a = 1, L — 100. 



3. Effect of a sampling cut off. 

We would like now to point out a very simple way to violate the scaling form (|3^), by the only numerical procedure 
traditionally used in the computation of (n) . This effect is closely related to the anomalous scaling of the Lee- Yang 
zeroes since it appears for a critical exponent r > 1. 
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In numerical simulations, one computes an empirical distribution P'^^ [n^uj) — — j!/^ where Af is the total number 
of avalanche observed during the simulation and M{n, oj) is the number of times an avalanche of size n was observed. 
This number is a random variable depending, say, on the initial condition(s), or more generally, on the seed w used in 
the random generator. However, one expects the system to be ergodic in a strong sense such that P^^^{n, uj) — > PL{n) 
as TV ^ cxo for generic choices of lo (see for details.). However, since M is finite, there exists wild fluctuations 
in the tail of the distribution. Furthermore, the avalanches such that Phin) < jj- have a small (though non zero []) 
probability to be observed in a numerical simulation and this probability decreases as n increases. Obviously, there 
are several methods such as smoothing or binning, allowing to reduce the effects of noise in the tail. 

There exists however another, more subtle effect. In all the examples of numerical computations that we have 
found in the SOC literature, the value of Af is fixed, independently of the system size. This induces a pathological bias 
affecting the extrapolation to the thermodynamic limit whatever the method used to analyze the empirical distribu- 
tion. In particular a violation of finite-size scaling can be observed on P'^^ {n,Lo) even if the theoretical probability 
Phin) obeys (^^. When TV is kept fixed while L increases, the estimation of the maximal value that the random 
variable can take (defining the exponent /3) is more and more biased. Indeed, while the true diverges as i ^ oo, 
the empirical value ^|f^(w) converges to a constant. Consequently, the probability distribution extrapolated to the 
thermodynamic limit from the empirical distribution is biased. The aim of this section is to analyze this effect, not 
discussed in the literature. 



. Call F^in) = J2k=i PL{k) the corresponding 
repartition function. Assume now that we perform a finite sampling of the probability distribution with TV trials 
Xi, . . . Xj^ where X^, i = 1 . . .M are independent |^, identically distributed random variable, with probability P^in). 
Assume furthermore that TV is fixed independently of L. Call ^2^^ ~ maxj^fe, k — \ . . .TV} the maximal value 
observed in the finite sampling. The repartition function of the random variable Cl^a/" 

Ff^{x). Its average is given 

by: 

E[W=iL-Y.PL{n) (47) 

ri=l 

Clearly, were L fixed while TV — ^ oo, then would £'[^2^^] converges to ^l- In fact the ergodic theorem gives a 
stronger statement, namely f^^j^ ~* almost-surely in this case. This essentially means that for sufficiently small 
L, E[^'^^^] gives a good estimate of ^l. On the other hand, if TV is fixed while L increases, the correction term 
"Y^^n^i Fif {^) in (0) becomes more and more important leading to a wrong estimation of ^l. To be more precise fix 
a value y e]0, 1[, such that is considered as non negligible as soon as F^{n) > y. Hence y is somehow arbitrary 
here (say close to 0). Since F^{n) is strictly increasing the equation F^{x) — y <^ Fl{x) — y^ has a unique solution 



'^However, since the largest non zero value of P'^^{n,uj) is jj, the events such that Phin) « jj', even when observed, are 
given an incorrect probability by the numerical procedure. The discrepancy with the theoretical value increases as n increases. 

^We assume here that the trials are independent for simplicity. In SOC models, the avalanches are not independent though, 
for finite L, the correlation decay can be fast (it is exponential in the finite size Zhang model). 



Assume therefore that Phin) is like (E9) with a cut off obeying (30 
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XL = XL^y), \/y £ ]0, 1[. If L is small (or if TV is sufficiently large) xl ^ and therefore F^{n) is essentially non 
zero for n ^ S^^. In this case the term X]n=i -^ifi''^) in ( p7| ) is negligible compared to f^. 

On the other hand, xl is bounded from above by a finite value x such that F*{x) = -ft'X]n=i '^"^ = 2/"^' where F* 
is the repartition function of the limiting probability P*{n). Consequently, as L increases, for any y g]0, 1[, ^^-^^^ — > 0. 
Hence, the function F^{x) is non negligible on a larger and larger interval, whose length scales like ^l- Therefore, the 
sum 

E^ti FlH in @ becomes more and more important, yielding a decrease in the expectation of the empirical 
maximum. 

In the range of L values where this effect starts to manifest one has xl ^l- Hence, Y^^=i Y^^=xl ^ifi''^) 

and Ffj' {n) has only to be estimated in the interval [a;^,^^]. Furthermore, Fi^{n) = 1 — X]i=n-^i('^) 1 — 
Cl /i"" u-^''gi-^)du. When xl is close to ^l, /j"' u"^ g{f^)du ~ (^^ - nj^l'' g{a) for n e [xL,iL\- Furthermore, 
in this range, 1 — Cl^^^{^l — n)9{o) is small compared to 1. Hence, F^ {n) ~ 1 — MClS.^^ — n)g{a). The equa- 
tion F^{x) = y has therefore an approximate solution xl ^ £,l ~ j\fc~[g[a) ■ Then J2n=i ^ifi^) Yl^n=XL ^'^) 
can be roughly approximated by a linear interpolation giving Yl^n=XL ^ifi''^) ^ ^^^-i^ ~ N'cT^ta) ^l- 



Consequently, the empirical expectation is, in this approximation: 



Since r > 1 the correction term increases as L grows. It eventually becomes of the same order as ^l, but when 
L increases one has to add higher order corrections to eq. (|8|). On the other hand, were t < 1, then would the 
correction term become negligible as L grows. 

The L value where the effect starts can be estimated by : 



Consequently, this effect is more prominent when /3(t — 1) is larger. 

The ratio aL = — i^ therefore not equal to a constant a as it must be, but is L dependent (see Fig. (o 
a)). Clearly, for sufficiently large L the corresponding probability violates the finite size scaling and the data collapse. 
Furthermore, the scaling of the moments is also affected by this effect. Indeed the moments are obtained empirically 

ex i"^'' ex 

from the formula m']^^(q) — X]n=i P^^^ (n) and the scaling exponents are extrapolated from the formula: 

^-P(^) ^ Hm M^S^ 

L^oo log(i) 

When L is sufficiently small, ^^^^ ^ aL^^ since the correction due to the finite sampling has essentially no effect. 
Then one obtains the right scaling exponent a{q) from the data. However, when L increases, one observes a spurious 
deviation of the curve mL{q) from the theoretical value. This effect is illustrated Fig. ^ in the case r ~ 1.9, /3 ~ 2,j — 
2.5, a = 2, g{x) — e^^^ where two samples with TV — 10® and TV — 10* were generated. Note that the effect is more 
prominent when q increases. 
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This computation shows therefore that for the values r > 1 numerical problems appears, induced by the finite size 
sampling. Not only the fluctuations of Clj^ increase but also the averaged value is biased. To our opinion, 

the estimation of the correct is the main problem in analyzing the data from SOC simulations. 

The analysis of the Lee- Yang zeroes for relatively small sizes can however give a fairly good estimate of the values 



a, [3 allowing to extrapolate to larger size. Indeed, eq. (40) suggests that the argument oitL{k) is not too sensitive 

''"^P-^ ). Hence, it 



to the fluctuations of (compared to the fluctuations of the moments which are of order ^Cl^a/") 
is possible to find the values of a, We give an example Fig.^ where the empirical data are the same as those used 
for the computation of the moments. As for the truncated power law, we found a slight deviation for the first zero. 
Interpolating the values from fc = 1 . . . 10 we found a = 2.094 ± 0.006 ± 0.09, /3 = 2.002 ± 0.003 which gives quite a 
good estimate. 
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FIG. 7. a) Argument of the first Lee- Yang zeroes of an empirical distribution generated from TV = 10* samples, b) Values 
of a,/3 extrapolated from a). 



The real part of Lee- Yang zeroes and consequently, the argument, is more sensitive to non extensive sampling effect. 
An analytical expression is obtained if one modifies the equations (Ea,E3^,E4h by replacing a by aL- The argument of 



20 



t^^k) writes now: 



Arg{tL{k)) ^ arctan 



2fc7r 



log(£-I|l^) + 1 log (log2(i-£|(^) + 4P7r2) 



(50) 



The finite sample effect is illustrated Fig. || for the first zero. One observes a deviation from the real curve for 
TV — 10^. This can be used as an empirical way to define the L where the empirical distribution is not biased. 
Note that the curve of the argument of the first empirical zero obtained for TV = 10^ follows the theoretical curve 
with a good accuracy. This shows that the determination of the zeroes is robust with respect to fluctuations in the 
coeflticients of the polynom (^. This is in fact an easy consequence of the implicit function theorem. 



0.35 



Theoretical value. 

Empirical value: /V=1 O^. 

& = Empirical value: /V=10^. 




20 30 40 50 60 70 80 90 1 0O 110 120 130 1 40 1 5Q 



FIG. 8. Effect of a sample cut off on the argument of the zeroes in the t plane. 



Looking at fig. (|p) and fig. (||) one could argue that the moments are less sensitive than the Lee- Yang zeroes to a 
size independent sampling. However, the sensitivity of the moments increases with the degree q. Therefore, to detect 
this effect, one has to compute the moments for high q. This can clearly cause numerical problems. On the other 
hand, the Lee- Yang zeroes integrate informations about all order moments and the sensitivity can be detected easily. 
We believe that this kind of analysis is prior to any investigation concerning in particular the multifractal nature of 
the scaling. 

Our conclusion is therefore twofold. Firstly, some cautions are required when increasing L to extrapolate the 
thermodynamic limit. If J\f is kept fixed, too large L will give wrong estimations. In this case at least, the bigger is 
not the best. This effect can be compared to the critical slowing down in the literature about critical phenomena. 
However, to the best of our knowledge we don't know any example in the SOC literature where this effect has been 
discussed. Secondly, the Lee- Yang zeroes give nevertheless useful informations. They can be used to determine the 
range of L values where the data are not too affected, and, in this range, the simple scaling of the imaginary part 
allows to determine the scaling of ^l- This can used with other methods such as binning, or smoothing to determine 
the exact distribution from the empirical one. 
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4- Other scaling. 



In this section we investigate briefly the other scaHng forms discussed in section II B . Our main conclusion is that 
the Lee- Yang zeroes are highly sensitive to the changes in the scaling form. This remark opens perspective to develop 
a general theory allowing to extrapolate the characteristic of probability distributions from the behavior of the Lee- 
Yang zeroes of the empirical generating function. However, we have not yet been able to provide an equivalent of 
the analytic forms (|2^^) that would be helpful to properly extract the features of the probability distribution from 
the Lee- Yang zeroes. The development of such a general framework is under investigation and will be published in a 
separated paper. We give a few numerical examples flg. Oa, b. 



We investigated first the case ( |33| ) where Fl{x) — Y{x — Pl)- In this case, the computations done in section [lA 
essentially hold, with a /3 depending on L. In particular, eq. ( p^ ) suggests that the L dependence of /3 should be 
detected on the argument of ZL{k), 9L{k). In fig. ^, we plotted 0l{5) in the case Pl = /?(! — x)' where /? = 2, r = 1.9, 
and L = 20 . . . 120. The theoretical prediction is 9l{5) = -^j^z^- We tried a fit of the form f{x) = -jj^^^r^, where 
a, (3 are the fit parameters. We found a = 1.014(8) ± 0.0006,/? = 2.013(9) ± 0.0001, which is quite satisfactory. 

We also tried a more general form for Fl{x) with a non linear Fl converging to a step function as L — > oo: 

FL{x) = -{l + iB.nh{L^^[x-P-^^y- x<P-^^ (51) 

In our simulations, /3 = 2, r = 1.9, ai — 0.1, a2 = 1. controls the rate of approach of Fl{x) to the step function 
as L grows. A small ai gives a slow convergence and therefore has an effect up to very large L. 

The result for the argument 6'l(5) is also represented fig. ||a. We note that the curve is indistinguishable from the 
previous case and therefore the imaginary part of is not sensitive to the non linear effect of the tanh, but gives the 
right a2- On the other hand, we noted that the argument of is sensitive to the non linear effect (fig. Sb). 



For the multifractal case we studied the case where the multifractal spectrum has the form f{x) = C — tx — ax^. 
This is the lowest degree non linear form of / compatible with (i),(ii) and with the convexity of /. The values of 
a — I, P — 2,T = 1.9 were the same as for the previous examples. We observe (Fig. ^) that f?L(fc) is not sensitive to 
the multifractality and gives therefore the right a, /3. Consequently, our method to estimate the degree is also valid 
for a multifractal distribution (^). On the other hand Arg{tL{k)) is modified for a multifractal distribution, but we 
are not yet able to analyze this variation. 

All the results are depicted fig. ||a (6Il(5) = Im{tL{5))) and fig. |b (Arg(tL(5))). 
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FIG. 9. Fig. Argument 9^(5) for the truncated power law (POW); a truncated power law with L dependent (3 (V); a 
probability distribution of the form ( ^3|) with Fl{x) given by (NL) and a multifractal distribution (MF).Fig. |^. Argument 

of fL(5). 



III. AN EXAMPLE OF SOC MODEL. 

In this section, we study an example of Lee- Yang zeroes computation for a SOC model, the Zhang model, defined 
as follows. 

Let A be a d-dimensional sub-lattice in Z'', taken as a square of edge length L for simplicity. Call N = ^A = L'*, 
and let ^A be the boundary of A, namely the set of points in Z'' \ A at distance 1 from A. Each site i e A is 
characterized by its "energy" Xi, which is a non-negative real number. Call X = {Xi}.^^ a configuration of energies. 
Let Ec be a real, strictly positive number, called the critical energy, and A4 = [0,Ec[^ . A configuration X is called 
"stable" if X G and "unstable" otherwise. If X is stable then one chooses a site i at random with probability 
and add to it energy ^(excitation). If a site i is over critical or active (X^ > Ec), it loses a part of its energy in equal 
parts to its 2d neighbors (relaxation). Namely, we fix a parameter e G [0, 1[ such that the remaining energy of i is 
eXi, after relaxation of the site i, while the 2d neighbors receive the energy '■"'""'^j"^' . Note therefore that the energy 
is locally conserved. If several nodes are simultaneously active, the local distribution rules are additively superposed, 
i.e. the time evolution of the system is synchronous. The succession of updating leading an unstable configuration to 
a stable one is called an avalanche. The energy is dissipated at the boundaries of the system, namely the sites of ^A 
have always zero energy. As a result, all avalanches are finite. Consequently, whatever the observable n, < oo for 
finite L. The addition of energy is adiabatic. When an avalanche occurs, one waits until it stops before adding a new 
energy quantum. Further excitations eventually generate a new avalanche, but, because of the adiabatic rule, each 
new avalanche starts from only one active site. It is conjectured that a critical state is reached in the thermodynamic 
limit. 

Though it has long been believed that the Zhang model obeys finite size scaling (pO|), a recent paper revised this 
point of view and claimed that the Zhang model does not even have a multifractal scaling (but no alternative scaling 
was proposed [|l5|). We will not solve this debate in this paper. Rather we will come to two conclusions. Firstly, 
because of high sensitivity of the model to the sample cut off (fig. O) , one has somehow to relativise the conclusions 



23 



about the scaling obtained from the numerical simulations. This shows that to draw any reliable conclusion on the 
scaling one has to increase the sample with the system size (e.g. like ). This will clearly be rapidly intractable even 
for the fastest computers. Secondly, the Lee- Yang zeroes gives rather reliable extrapolations provided the size L is 
not too large. 



We computed the empirical probability distribution of avalanche sizes P^^{s) where the size is the total number 
of relaxing sites during one avalanche. We did our simulations for lattice sizes from L = 10 to L = 55 in two 
dimensions, with Ec = 2.2, e = 0.1 with a statistics over J\f = 10^ and TV = 10* avalanches. Consequently, J\f was 
fixed independently of L as usually done in SOC numerical simulations. 

We first present Fig. [lo| a the experimental zeroes in the z complex plane for L = 30, A/" = 10^. In order to see 
the effect of the noise, we also computed the zeroes of a smoothed version of the empirical probability distribution 
(fig. (p^)). The smoothing method uses a binning procedure, followed by a spline extrapolation, allowing to fill 
the "holes" existing in the empirical probability distribution. These holes corresponds to events that didn't happen 
during the trial and consequently are given a zero probability. In the numerical computation of the zeroes, these holes 
corresponds therefore to vanishing coefficients in the polynom , that produce problems in the convergence of most 
of the root finding algorithms. Our numerical procedure seems however to be robust with respect to this effect. 
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FIG. 10. a) Zeroes in tlie z plane for L = 30, e = 0.1, TV = 10^, Ec = 2.2.. The zeroes of the experimental, TV = 10^, and 
smoothed probability distribution are represented, b) Examples of empirical and smoothed probability distributions used in 
the computation of the zeroes 



As argued all along in this paper, the argument OL{k) provides a way to determine the exponent (3 characterizing 
the maximal avalanche size. We plot fig. |ll|a the 6Il(A:)'s for fc = 1 . . . 5 versus L. We note that dL{k) is quite robust 
to noise and gives therefore a reliable way to measure /3, a. We used a fit form ^j^, where a, (3 are fit parameters. We 
found a slight k dependence for the first zeroes (as expected from eq. ( p3| ) if one assumes FSS). We plotted fig. [l]b 
the extrapolated a, (3 versus k. For fc > 3 these value seem to stabilize around a = 0.62 ± 0.07 and /3 — 2.59 ± 0.04. In 
the finite size scaling ansatz (3 and r are related by (3{2 — t) = 2 [ p5[ and r — 1.253 is known from the renormalization 
group analysis. Therefore the predicted value for (3 is 2.667. Despite the smallness of the L we considered, the 
computed value is not too far from the predicted one. However, an accurate determination of j3, r and a precise check 
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of FSS would demand somewhat larger size systems, that we were unable to generate for this illustration. (Note that 
the main problem is not the computation of the zeroes since there exist quite fast and precise root finding algorithms, 
but the generation of P^^^ itself.). 




10 k=1 2 4 6 8 10 

FIG. 11. a. Argument of the 5 first Lee- Yang zero in the z plane versus L, for Ec — 2.2, e = 0.1. The full lines corresponds 
to a fit. b. Plot of the fit parameters a,/3 versus k . 



Finally, we investigate the scaling of the argument tL{k) and a possible size independent sampling effect. The main 
difficulty here are the wild fluctuations in ^'j^^ . Indeed, the real part of t^^k) is more sensitive than the imaginary part 
to these fluctuations and consequently Arg{tL(k)) as wild fluctuations. Only the first zeroes seem to be robust to this 



exp{L) 



effect. As discussed in section [IB 3 these fluctuations in are intrinsic to the empirical computation of 
and cause problems in the extrapolation to the thermodynamics limit, whatever the method used. We plotted Fig.[l^ 
the argument of the 2 first Lee- Yang zeroes in the t plane, versus L. Our simulations suggests that the Zhang model 
is sensitive to the size independent sampling. Note in particular that the values of L that we used in our computation 
are quite smaller that the ones found in the literature and the effect is already significant for JV < 10®. 
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FIG. 12. Argument of Lee- Yang zeroes tL(l), tL(l) for Ec = 2.2, e — 0.1. In color full lines are drawn the interpolation of 
the empirical curves obtained from the sampling with Af = 10*. 



This shows clearly that a re estimation of the conclusions drawn from the numerics in the Zhang model (and also, 
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maybe, for some other SOC models) should be done in light of this effect. On the other hand, our approach suggests 
that there may be no need to go to gigantic sizes provided the finite lattice size effects are carefully handled. From 
this point of view, analytic formula like (Q) might provide a way to analyze these effects. 

IV. CONCLUSION. 

In this paper, we have shown that the finite size study of the SOC like probability distributions leads to similar Lee 
Yang or Fisher phenomenon as in statistical physics models near the critical point. This implies that the convergence 
of the SOC state to a critical state with power law statistics can be analyzed in a similar way as equilibrium statistical 
mechanics. More precisely, the way the zeros of the partition function accumulate on the real axis, when the size of 
the system grows up, provides relevant informations on the critical structure of the observed system. In particular, it 
permits to measure useful critical indexes of the underlying theory. 

Moreover, we have shown that the size of the SOC models power exponent, r > 1, leads to a comprehensive violation 
of the standard scaling laws. We give a approximate theory of this effect well confirmed by numerical simulations. It 
is the same characteristic which leads to a specific sensibility of the SOC numerical experiments to size independent 
sampling effects. We studied carefully this effect on extrapolation to the limit L ^ oo and show that it could possibly 
mimics important effects such as multifractality. We notice that the argument of the zeros in the t = log{z) plane is 
a good test of this effect. 

On one other hand, we show that the arguments of the first zeros in the z plane of the generating function, G{z) is 
rather insensitive to these effects, statistically robust, and provide a nice way to compute the SOC a and /? parameters. 
Using the standard Kadanoff et al. scaling form [ p!3[ , we verify that the parameter's values as extracted from numerical 
simulations where in good agreement with the theoretical input of the model. This last result gives us some confidence 
to extract the values of these parameters from Zhang's model numerical data. Notice that these results have been 
extracted from medium range simulations. This shows up once more the power of the finite size analysis of the critical 
phenomenon. 

This paper is (with |^) a first step toward a scaling theory of SOC system from the behavior of the Lee- Yang 
zeroes. The next step would be the definition of the exponents characterizing the approach to criticality, like the 
exponents a, (3, 7 in statistical mechanics and their link to the scaling of the zeroes. 
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